
(*-------------------------------------------------------------------------------------------------*)
(*                                                                                                 *)
(*                           M A T H E M A T I C A   -   S C R I P T                               *)
(*                                                                                                 *)
(*           zum Berechnen aller Ableitungen der 1D Zylindergeoemetrie RHD-Gleichungen             *)
(*                                                                                                 *)
(*                                                                                                 *)
(*     Benutzung:                                                                                  *)
(*     - gewuenschte Gleichung ein-kommentieren                                                    *)
(*     - mit  "math"  Mathematica auf der command-line starten                                     *)
(*     - vorliegende Datei mittels   "<< 1DD-RHD.m"  evaluieren                                    *)
(*     - FORTRAN-Code aus output.f90 in die entsprechende DG_... Datei kopieren                    *)
(*                                                                                                 *)
(*                                                                                                 *)
(*-------------------------------------------------------------------------------------------------*)



(*-------------------------------------------------------------------------------------------------*)
(*                           E i n l e s e n   d e r   P a c k a g e s                             *)
(*-------------------------------------------------------------------------------------------------*)

    (* Setzen des Arbeitsverzeichnisses zum Einlesen der Packages *)

  SetDirectory["./packages"];

    (* Beim Einlesen der Packages Variablenhirachie beachten ! *)

  << primVar.m;         (* Definition der primaeren Variablen *)
  << zz.m;              (* Definition des zeitlich zentrierten Variablensatzes *)
  << zint.m;            (* Definition der Variablen fuer die Integration in z-Richtung *)
  << matfct.m;          (* Definition der Materialfunktionen als Funktionen *)
  << advec.m;           (* Definition der advektierten Variablen als Funktionen *)

  << geomet.m;          (* Definition der geometrischen Hilfsvariablen - Volumina & Fluesse *)
  << grid.m;            (* Definition der Hilfsvariablen fuer die Gittergleichung *)
  << visco.m;           (* Definition der Hilfsvariablen fuer die Viskositaet *)

  << sekVar.m;          (* Ersetzungsfunktionen fuer Hilfsvariable *)
  << ForVar.m;          (* Ersetzungen der mathamatica Variablennamen durch FORTRAN Variablennamen *)


    (* Zuruecksetzen des Arbeitsverzeichnisses *)

  ResetDirectory[ ];


    (* weitere vorbereitende Kleinigkeiten... *)

  pagewidth = 130;      (* Zeilenlaenge im output-file *)
  Off[];                (* Unterduecken unnoetiger Fehlermeldungen *)
  Run["rm output.f90"]; (* Loeschen des alten output-files *)



  Print["Ableitungen der ..."];
(*-------------------------------------------------------------------------------------------------------*)
(*                       D e f i n i t i o n    d e r    G l e i c h u n g e n :                         *)
(*-------------------------------------------------------------------------------------------------------*)
(*
  Print[">>> Gittergleichung"];
  Equation = "MR";

  G = Xntmpi/XResi - Xntmpim1/XResim1;
  *)
(*
  Print[">>> Gittergleichung Randbedingung"];
  Equation = "MR";

  G = Xntmpi - Xntmpim1;
  *)
(*
  Print[">>> Gittergleichung Lagrange"];
  Equation = "MR";

  G = XSfluxi;
  *)

(*-------------------------------------------------------------------------------------------------------*)
(*
  Print[">>> Kontinuitaetsgleichung"];
  Equation = "MD";

  G = XSvoli rhoi - SvolAi rhoAi \
    + XSfluxip1 Fktrhoadvip1 - XSfluxi Fktrhoadvi;
  *)

(*-------------------------------------------------------------------------------------------------------*)
(*
  Print[">>> Energiegleichung"];
  Equation = "ME";

  G = XSvoli ei rhoi - SvolAi eAi rhoAi \
    + XSfluxip1 Fkteadvip1 - XSfluxi Fkteadvi \
    + sqrt2pi FktHpZi FktPgas0Zi z2 pi (XRZip1 XurZip1 - XRZi XurZi) tst \
    - z4 pi FktOPAplkZi XrhoZi (XJZi - StBolz/pi FktTgasZi^4) XSvolZi tst \
    - z12 (XmuQri XUtermri^2 + XmuQpenei XUtermpi^2) XSvolZi tst \
    - Lstar firr HpZflagi ( (FktHpZip1-Hstar)/XRZip1 - (FktHpZi-Hstar)/XRZi ) tst \
    + z8 pi muz XJZi XSvolZi tst;
  *)

(*-------------------------------------------------------------------------------------------------------*)
(*
  Print[">>> Bewegungsgleichung radial"];
  Equation = "MUr";

  G = XVvoli uri rhoi - VvolAi urAi rhoAi \
    + XVfluxi Fktmomradvi - XVfluxim1 Fktmomradvim1 \
    - z12 (XupZi^2 XSvolZi + XupZim1^2 XSvolZim1) / XRZi XrhoZi tst \
    + z2 pi XRZi sqrt2pi (FktHpZi FktPgas0Zi - FktHpZim1 FktPgas0Zim1) tst \
    + grav Mstar / XRZi^2 XrhoZi XVvolZi tst \
    - z2 pi / clight XHZi (FktOPArosZi XrhoZi XSvolZi + FktOPArosZim1 XrhoZim1 XSvolZim1) tst \
    - z12 pi / XRZi ( XmuQri (XRZi^2 + XRZip1^2) XUtermri - XmuQrim1 (XRZi^2 + XRZim1^2) XUtermrim1 ) tst;
  *)

(*-------------------------------------------------------------------------------------------------------*)
(*
  Print[">>> Bewegungsgleichung phi"];
  Equation = "MUp";

 G = XSvoli upi rhoi z12 (Ri + Rip1) - XSvolAi upAi rhoAi z12 (RAi + RAip1) \
   + XSfluxip1 XRZip1 Fktmompadvip1 - XSfluxi XRZi Fktmompadvi \
   - z12 pi ( XmuQpip1 XRZip1^2 XUtermpip1 - XmuQpi XRZi^2 XUtermpi ) tst;

  *)

(*-------------------------------------------------------------------------------------------------------*)
(*
  Print[">>> Strahlungsenergiegleichung"];
  Equation = "MJ";

  G = XSvoli Ji - SvolAi JAi \
    + XSfluxip1 FktJadvip1 - XSfluxi FktJadvi \
    + clight z2 pi (Rip1 Hip1 - Ri Hi) tst \
    + feddi Ji z2 pi (Rip1 urip1 - Ri uri) tst \
    - z12 Ji (z3 feddi - z1) XSvoli tst (Ri uri + Rip1 urip1) / (Ri^2 + Rip1^2) \
    + clight FktOPAplki Xrho0i (Ji - StBolz/pi FktTgasi^4) XSvoli tst;
  *)

(*
  Print[">>> Strahlungsenergiegleichung: J=S"];
  Equation = "MJ";

  G = Ji - StBolz/pi FktTgasi^4;
  *)

(*-------------------------------------------------------------------------------------------------------*)
(*
  Print[">>> Strahlungsflussgleichung"];
  Equation = "MH";

  G = XVvoli Hi - VvolAi HAi \
    + XVfluxi FktHadvi - XVfluxim1 FktHadvim1 \
    + clight z2 pi Ri (feddi Ji - feddim1 Jim1) tst \
    + clight z14 ( Ji (z3 feddi - z1) XSvoli + Jim1 (z3 feddim1 - z1) XSvolim1 ) / Ri tst \
    + Hi pi Ri (urip1 - urim1) tst \
    + clight z12 (FktOPArosi Xrho0i XSvoli + FktOPArosim1 Xrho0im1 XSvolim1) Hi tst;
  *)

(*
  Print[">>> Abstrahlrandbedingung"];
  Equation = "MH";

  G = Hi - z12 (XSvoli Ji + XSvolim1 Jim1) / XVvoli mur;
  *)

(*-------------------------------------------------------------------------------------------------------*)
(*                                          T E S T  -  A r e a                                          *)
(*-------------------------------------------------------------------------------------------------------*)
(**)
  Print[">>> TEST-Area"];
  Equation = "MT";

  G =  Ri + rhoi + uri + upi + ei + Ji + Hi;
(**)

(*-------------------------------------------------------------------------------------------------------*)
(*              B e r e c h n u n g   u n d   A u s g a b e   d e r   A b l e i t u n g e n              *)
(*-------------------------------------------------------------------------------------------------------*)


    (* Aufstellen der Jacobi-Matrix & Anwednung der Ersetzungregeln *)

  DG = Table[baksubzz[baksubzint[baksubgeomet[baksubgrid[baksubvisco[D[G, primvar[[(i-1)*np + j]]]]]]]], {i,1,nv}, {j,1,np}];
  TableForm[DG];


    (* Ausgabe der Ableitungen in output.f90 *)

  Do[out = OpenAppend["output.f90", PageWidth -> pagewidth];
      If[ToString[DG[[ii, jj]]] != "0", \
         Write[out, "DG(",Equation,",",VarMarker[ii],",",DGidxMarker[jj],",i)","=", \
                    fortranvariable[ToString[FortranForm[zkonstante[N[DG[[ii, jj]] ]]]]] \
              ] \
        ];
     Close["output.f90"], {ii, 1, nv}, {jj, 1, np} \
    ]


(*                    fortranvariable[ToString[FortranForm[zkonstante[Simplify[N[DG[[ii, jj]] ],TimeConstraint->1]]]]] \   *)

(*-------------------------------------------------------------------------------------------------------*)
(*                                 P o s t  --  P r o c e s s i n g                                      *)
(*-------------------------------------------------------------------------------------------------------*)

  Run["./postprocess"];




